An exploration of effective and useful visualisation of optical satellite images

Make your satellite image pop!

In [1]:
import re

import numpy as np
import matplotlib.pyplot as plt

import boto3
import rasterio
from rasterio.plot import show, show_hist
from skimage.exposure import equalize_adapthist
In [3]:
session = boto3.Session(profile_name='default')
In [4]:
scenes = [
    'LC08_L1TP_205021_20180625_20180704_01_T1',
    'LC08_L1TP_205021_20180727_20180731_01_T1'
]
In [5]:
def get_aws_band(scene, band, profile=False):
    path, row = re.match(r'LC08_L1TP_(\d{3})(\d{3})', scene).groups()
    prefix = f'c1/L8/{path}/{row}/{scene}'
    
    with rasterio.Env(session=session):
        with rasterio.open(
            f's3://landsat-pds/{prefix}/{scene}_B{band}.TIF') as src:
            
            if profile:
                return src.profile
            
            return src.read(out_shape=(src.height//10, src.width//10))
In [6]:
stacked_scenes = []
for scene in scenes:
    red = get_aws_band(scene, 4)[0]
    green = get_aws_band(scene, 3)[0]
    blue = get_aws_band(scene, 2)[0]

    stacked_scenes.append(np.stack([red, green, blue]))

Let's plot our downloaded image

We will use rasterio's plotting function (which wraps matplotlib)

In [11]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(stacked_scenes[0], title="Scene 1", ax=axes[0])
show(stacked_scenes[1], title="Scene 2", ax=axes[1])
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).
Out[11]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa55a42aeb8>
In [12]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,65535)
axes[1].set_xlim(0,65535)
axes[0].set_ylim(0,140000)
axes[1].set_ylim(0,140000)
show_hist(
    stacked_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    stacked_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])
In [13]:
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif
Driver: GTiff/GeoTIFF
Files: data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif
Size is 8061, 8151
Coordinate System is:
PROJCS["WGS 84 / UTM zone 30N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32630"]]
Origin = (341085.000000000000000,6318015.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  341085.000, 6318015.000) (  5d36'53.37"W, 56d58'41.99"N)
Lower Left  (  341085.000, 6073485.000) (  5d28'16.42"W, 54d47' 0.10"N)
Upper Right (  582915.000, 6318015.000) (  1d38' 6.39"W, 56d59'53.56"N)
Lower Right (  582915.000, 6073485.000) (  1d42'36.48"W, 54d48' 6.03"N)
Center      (  462000.000, 6195750.000) (  3d36'28.13"W, 55d54'20.48"N)
Band 1 Block=8061x1 Type=UInt16, ColorInterp=Gray
Band 2 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 3 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 4 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 5 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 6 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 7 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 8 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 9 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 10 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 11 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 12 Block=8061x1 Type=UInt16, ColorInterp=Undefined
In [14]:
!gdalinfo data/LC08_L1TP_205021_20180625_20180704_01_T1/full_stack.tif | grep "Band"
Band 1 Block=8061x1 Type=UInt16, ColorInterp=Gray
Band 2 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 3 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 4 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 5 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 6 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 7 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 8 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 9 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 10 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 11 Block=8061x1 Type=UInt16, ColorInterp=Undefined
Band 12 Block=8061x1 Type=UInt16, ColorInterp=Undefined

Ok, let's just make it 8bit, that'll work, right...

Most plotting tools work with 8 bit images (we have 16 bit):

  • Red; 256 tonal levels
  • Green; 256 tonal levels
  • Blue; 256 tonal levels

In [15]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(stacked_scenes[0].astype(np.uint8), title="Scene 1", ax=axes[0])
show(stacked_scenes[1].astype(np.uint8), title="Scene 2", ax=axes[1])
Out[15]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa569ab7128>
In [16]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,50000)
axes[1].set_ylim(0,50000)
show_hist(
    stacked_scenes[0].astype(np.uint8), bins=50, lw=0.0, stacked=False,
    alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    stacked_scenes[1].astype(np.uint8), bins=50, lw=0.0, stacked=False,
    alpha=0.3, histtype='stepfilled', title="Histogram", ax=axes[1])

That's pretty psychedelic, but not helpful

We have to be a bit smarter than that, we need to rescale!

In [17]:
def linear_rescale(image, in_range=[0, 65535], out_range=[1, 255]):
    """
    Linear rescaling
    """

    image = np.clip(image, in_range[0], in_range[1]) - in_range[0]
    image = image / float(in_range[1] - in_range[0])

    return (
        image * (out_range[1] - out_range[0]) + out_range[0]
    ).astype(np.uint8)
In [18]:
linear_rescaled_images = [linear_rescale(stack) for stack in stacked_scenes]
In [19]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
Out[19]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa5539a90b8>
In [20]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
    linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])

OK, Landsat 8 doesn't actually make use of the full 16 bits, so we will rescale within the range we have reasonable values:

In [21]:
linear_rescaled_images = [
    linear_rescale(stack, in_range=[0, 14000]) for stack in stacked_scenes
]
In [22]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(linear_rescaled_images[1], title="Scene 2", ax=axes[1])
Out[22]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa5537ff438>
In [23]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
    linear_rescaled_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    linear_rescaled_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])

NASA

That's not very pretty

We can make better use of the colour range available, so lets stretch it.

In [24]:
def percentile_stretch(image, lower_perc=2, upper_perc=98):
    mask = np.all(image != 0, axis=0).astype(np.uint8) * 255
    image_out = np.zeros_like(image).astype(np.uint8)
    for band in range(image.shape[0]):
        input_range = np.percentile(
            image[band][mask > 0],
            (lower_perc, upper_perc)
        ).tolist()
        
        image_out[band] = np.where(
            mask,
            linear_rescale(image[band],
            in_range=input_range,
            out_range=[0, 255]),
            0
        )
    
    return image_out
In [25]:
percentile_stretched_images = [
    percentile_stretch(stack) for stack in stacked_scenes
]
In [26]:
fig, axes = plt.subplots(1, 2, figsize=(30,15))

show(percentile_stretched_images[0], title="Scene 1", ax=axes[0])
show(percentile_stretched_images[1], title="Scene 2", ax=axes[1])
Out[26]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa55394f828>
In [27]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
    percentile_stretched_images[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    percentile_stretched_images[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])

Give me more detail

That worked nicely with our clearish scene, but in the cloudy one, well the clouds look really interesting.
Contrast is lacking elsewhere...

Introducing Contrast Limited Adaptive histogram Equalization

In [28]:
def clahe_equalize(image, kernel_size=2000, clip_limit=0.03):
    image_out = np.zeros_like(image).astype(np.float64)
    for band in range(image.shape[0]):
        image_out[band] = equalize_adapthist(
            image[band], kernel_size=kernel_size, clip_limit=clip_limit
        )

    return image_out
In [29]:
raster_equalized_images = [
    clahe_equalize(stack, clip_limit=0.05) for stack in stacked_scenes
]
In [30]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(
    (raster_equalized_images[0]* 255).astype(np.uint8),
    title="Scene 1", ax=axes[0]
)
show(
    (raster_equalized_images[1]* 255).astype(np.uint8),
    title="Scene 2", ax=axes[1]
)
Out[30]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa5536db048>
In [31]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,80000)
axes[1].set_ylim(0,80000)
show_hist(
    (raster_equalized_images[0]* 255).astype(np.uint8),
    bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
    title="Histogram", ax=axes[0]
)
show_hist(
    (raster_equalized_images[1]* 255).astype(np.uint8),
    bins=50, lw=0.0, stacked=False, alpha=0.3, histtype='stepfilled',
    title="Histogram", ax=axes[1]
)

We could also remove the atmosphere...kind of...

The scenes still look very different, so not great for applying to all scenes, let's think about this in a physical sense.

We could try and remove the atmosphere, and ideally with an image based method (simpler, quicker, and no need for additional external data)

In [32]:
def find_dark_object_value(arr):
    """Find the value of the dark object
    ie the first dark value that is not an outlier
    """
    preval = None
    step = arr.max() / 255.0
    for val in np.unique(arr)[:100]:
        if val == 0:
            continue
        if preval is not None and (val - preval) < step:
            break
        else:
            preval = val
    return preval


def dos(image):
    image_out = np.zeros_like(image).astype(np.uint8)
    
    # get the dark object for each band
    for i, band in enumerate(image): 
        dark = find_dark_object_value(band)
        image_out[i] = band - dark.astype(np.float)

    image_out[image_out < 0] = 0
    
    return image_out.astype(np.uint8)
In [33]:
prepped_scenes = [
    linear_rescale(stack, [0, 14000]) for stack in stacked_scenes
]
dos_scenes = [dos(prepped_scene) for prepped_scene in prepped_scenes]
In [34]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(linear_rescaled_images[0], title="Scene 1", ax=axes[0])
show(dos_scenes[0], title="Scene 1 - DOS", ax=axes[1])
Out[34]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa55353c940>
In [35]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(linear_rescaled_images[1], title="Scene 2", ax=axes[0])
show(dos_scenes[1], title="Scene 2 - DOS", ax=axes[1])
Out[35]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa568e31dd8>
In [36]:
fig, axes = plt.subplots(1,2, figsize=(30,15))

show(dos_scenes[0], title="Scene 1 - DOS", ax=axes[0])
show(dos_scenes[1], title="Scene 2 - DOS", ax=axes[1])
Out[36]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fa568d85160>
In [37]:
fig, axes = plt.subplots(1,2, figsize=(30,10))

axes[0].set_xlim(0,256)
axes[1].set_xlim(0,256)
axes[0].set_ylim(0,120000)
axes[1].set_ylim(0,120000)
show_hist(
    dos_scenes[0], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[0])
show_hist(
    dos_scenes[1], bins=50, lw=0.0, stacked=False, alpha=0.3,
    histtype='stepfilled', title="Histogram", ax=axes[1])

Export for display

In [38]:
def write_rendered_rgb(rgb, profile, outname_prefix):
    """
    Write a rendered RGB to GeoTIFF.
    """
    profile.update(dtype=rasterio.uint8, count=3, photometric='RGB')

    outname = f'data/{outname_prefix}_rendered.tif'

    with rasterio.open(outname, 'w', **profile) as dst:
        for k, arr in enumerate(rgb):
            dst.write_band(k+1, arr)
In [39]:
profile = get_aws_band(scenes[1], 4, True)
equalized_image_8bit = (raster_equalized_images[1] * 255).astype(np.uint8)
write_rendered_rgb(equalized_image_8bit, profile, "scene2")
In [41]:
!gdalinfo data/scene2_rendered.tif | grep Band
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue

Ta